*** Set working directory ***
cd "C:\..."

*********************** Load Data **********************************************
clear
graph drop _all
cap drop all

import excel "Data.xlsx", sheet("data_extended") firstrow

************************** Setup ***********************************************
* Choose impulse response horizon
local hmax = 36

* Time variable
gen monthly_date = mofd(date)
format monthly_date %tm
drop date
rename monthly_date date
sort date
tsset date
egen time = group(date)
drop if date == .

* Standardize MP Shock series - unit standard deviation
egen mps_std = sd(chg_sveny02)
gen mps = chg_sveny02/mps_std

sort date
tsset date

* Standardize GSCPI for 1998-2023 sample
rename gscpi gscpi_unscl
egen gscpi_mean = mean(gscpi_unscl)
gen gscpi_dm = gscpi_unscl - gscpi_mean
egen gscpi_std = sd(gscpi_dm)
gen gscpi = gscpi_dm/gscpi_std

* GSCPI mp interaction terms
gen gscpi_mps = L.gscpi*mps
gen gscpi_mps1 = L.gscpi*L.mps
gen gscpi_mps2 = L.gscpi*L2.mps

* Put variables in log form
gen cpi_inf = ln(cpi)
gen rgdp_growth = ln(ip)
rename unrt unrt_chg
rename retail retail2
gen retail = ln(retail2)

* SP500 in log form
gen stk = ln(sp)

* Endog variables
forvalues h = 0/`hmax' {
	gen stk_`h' = f`h'.stk
	gen ebp`h' = f`h'.ebp
	gen twotp`h' = f`h'.twotp
	gen tentp`h' = f`h'.tentp
}

************************** LP regressions **************************************
* S&P 500
eststo clear
cap drop b1 b2 u1a u2a d1a d2a u1 u2 d1 d2 Years1 Zero1
gen Years1 = _n-1 if _n<=`hmax'
gen Zero1 =  0    if _n<=`hmax'
gen b1=0
gen b2=0
gen u1a=0
gen u2a=0
gen d1a=0
gen d2a=0
gen u1=0
gen u2=0
gen d1=0
gen d2=0
gen stde1 = 0 

forv h = 0/`hmax' {
	
	newey stk_`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg l(0/2).retail l(1/2).stk, lag(`h')
replace b1 = _b[mps]                    if _n == `h'+1
replace u1a = _b[mps] + 1.645* _se[mps]  if _n == `h'+1
replace d1a = _b[mps] - 1.645* _se[mps]  if _n == `h'+1
replace u1 = _b[mps] + 1* _se[mps]  if _n == `h'+1
replace d1 = _b[mps] - 1* _se[mps]  if _n == `h'+1

replace b2 = _b[mps] + _b[gscpi_mps]                    if _n == `h'+1
matrix C1=(1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
matrix Vcomb1 =C1*e(V)*C1'         
mat SEcomb1 =cholesky(Vcomb1)      
scalar temp1 = SEcomb1[1,1]        
replace stde1 = temp1            

replace u2a = _b[mps] + _b[gscpi_mps] + 1.645* stde1  if _n == `h'+1
replace d2a = _b[mps] + _b[gscpi_mps] - 1.645* stde1  if _n == `h'+1
replace u2 = _b[mps] + _b[gscpi_mps] + 1* stde1  if _n == `h'+1
replace d2 = _b[mps] + _b[gscpi_mps] - 1* stde1  if _n == `h'+1

eststo
}


* EBP
cap drop b3 b4 u3a u4a d3a d4a u3 u4 d3 d4 Years2 Zero2
gen Years2 = _n-1 if _n<=`hmax'
gen Zero2 =  0    if _n<=`hmax'
gen b3=0
gen b4=0
gen u3a=0
gen u4a=0
gen d3a=0
gen d4a=0
gen u3=0
gen u4=0
gen d3=0
gen d4=0
gen stde2 = 0 

forv h = 0/`hmax' {
	
	newey ebp`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg  l(0/2).retail l(1/2).ebp, lag(`h')
replace b3 = _b[mps]                    if _n == `h'+1
replace u3a = _b[mps] + 1.645* _se[mps]  if _n == `h'+1
replace d3a = _b[mps] - 1.645* _se[mps]  if _n == `h'+1
replace u3 = _b[mps] + 1* _se[mps]  if _n == `h'+1
replace d3 = _b[mps] - 1* _se[mps]  if _n == `h'+1

replace b4 = _b[mps] + _b[gscpi_mps]                    if _n == `h'+1
matrix C2=(1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
matrix Vcomb2 =C2*e(V)*C2'         
mat SEcomb2 =cholesky(Vcomb2)      
scalar temp2 = SEcomb2[1,1]        
replace stde2 = temp2            

replace u4a = _b[mps] + _b[gscpi_mps] + 1.645* stde2  if _n == `h'+1
replace d4a = _b[mps] + _b[gscpi_mps] - 1.645* stde2  if _n == `h'+1
replace u4 = _b[mps] + _b[gscpi_mps] + 1* stde2  if _n == `h'+1
replace d4 = _b[mps] + _b[gscpi_mps] - 1* stde2  if _n == `h'+1

eststo
}


* 2 year TP
cap drop b5 b6 u5a u6a d5a d6a u5 u6 d5 d6 Years3 Zero3
gen Years3 = _n-1 if _n<=`hmax'
gen Zero3 =  0    if _n<=`hmax'
gen b5=0
gen b6=0
gen u5a=0
gen u6a=0
gen d5a=0
gen d6a=0
gen u5=0
gen u6=0
gen d5=0
gen d6=0
gen stde3 = 0 

forv h = 0/`hmax' {
	
	newey twotp`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg  l(0/2).retail l(1/2).twotp, lag(`h')
replace b5 = _b[mps]                    if _n == `h'+1
replace u5a = _b[mps] + 1.645* _se[mps]  if _n == `h'+1
replace d5a = _b[mps] - 1.645* _se[mps]  if _n == `h'+1
replace u5 = _b[mps] + 1* _se[mps]  if _n == `h'+1
replace d5 = _b[mps] - 1* _se[mps]  if _n == `h'+1

replace b6 = _b[mps] + _b[gscpi_mps]                    if _n == `h'+1
matrix C3=(1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
matrix Vcomb3 =C3*e(V)*C3'         
mat SEcomb3 =cholesky(Vcomb3)      
scalar temp3 = SEcomb3[1,1]        
replace stde3 = temp3            

replace u6a = _b[mps] + _b[gscpi_mps] + 1.645* stde3  if _n == `h'+1
replace d6a = _b[mps] + _b[gscpi_mps] - 1.645* stde3  if _n == `h'+1
replace u6 = _b[mps] + _b[gscpi_mps] + 1* stde3  if _n == `h'+1
replace d6 = _b[mps] + _b[gscpi_mps] - 1* stde3  if _n == `h'+1

eststo
}


* 10 year TP
cap drop b7 b8 u7a u8a d7a d8a u7 u8 d7 d8 Years4 Zero4
gen Years4 = _n-1 if _n<=`hmax'
gen Zero4 =  0    if _n<=`hmax'
gen b7=0
gen b8=0
gen u7a=0
gen u8a=0
gen d7a=0
gen d8a=0
gen u7=0
gen u8=0
gen d7=0
gen d8=0
gen stde4 = 0 

forv h = 0/`hmax' {
	
	newey  tentp`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg l(0/2).retail l(1/2).tentp, lag(`h')
replace b7 = _b[mps]                    if _n == `h'+1
replace u7a = _b[mps] + 1.645* _se[mps]  if _n == `h'+1
replace d7a = _b[mps] - 1.645* _se[mps]  if _n == `h'+1
replace u7 = _b[mps] + 1* _se[mps]  if _n == `h'+1
replace d7 = _b[mps] - 1* _se[mps]  if _n == `h'+1

replace b8 = _b[mps] + _b[gscpi_mps]                    if _n == `h'+1
matrix C4=(1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
matrix Vcomb4 =C4*e(V)*C4'         
mat SEcomb4 =cholesky(Vcomb4)      
scalar temp4 = SEcomb4[1,1]        
replace stde4 = temp4            

replace u8a = _b[mps] + _b[gscpi_mps] + 1.645* stde4  if _n == `h'+1
replace d8a = _b[mps] + _b[gscpi_mps] - 1.645* stde4  if _n == `h'+1
replace u8 = _b[mps] + _b[gscpi_mps] + 1* stde4  if _n == `h'+1
replace d8 = _b[mps] + _b[gscpi_mps] - 1* stde4  if _n == `h'+1

eststo
}


************************* Plot IRF's *******************************************
* S&P 500
twoway ///
(rarea u1a d1a  Years1,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u1 d1  Years1,  ///
fcolor(gs12) lcolor(white) lpattern(solid))  ///
(line b1 Years1, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero1 Years1, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_gdp1, replace

twoway ///
(rarea u2a d2a  Years1,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u2 d2  Years1,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b2 Years1, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero1 Years1, lcolor(black)), legend(off) ///
title("Tight Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_gdp2, replace

gr combine fig_gdp1 fig_gdp2, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("S&P 500", size(small)) ///

gr rename fig_ebp, replace


* EBP
twoway ///
(rarea u3a d3a  Years2,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u3 d3  Years2,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b3 Years2, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero2 Years2, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_cpi1, replace

twoway ///
(rarea u4a d4a  Years2,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u4 d4  Years2,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b4 Years2, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero2 Years2, lcolor(black)), legend(off) ///
title("Tight Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_cpi2, replace

gr combine fig_cpi1 fig_cpi2, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("EBP", size(small)) ///

gr rename fig_cp, replace


* 2 year TP 
twoway ///
(rarea u5a d5a  Years3,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u5 d5  Years3,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b5 Years3, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero3 Years3, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ur1, replace

twoway ///
(rarea u6a d6a  Years3,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u6 d6  Years3,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b6 Years3, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero3 Years3, lcolor(black)), legend(off) ///
title("Tight Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ur2, replace

gr combine fig_ur1 fig_ur2, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("Two-Year Term Premium", size(small)) ///

gr rename fig_aaa, replace


* 10 year TP 
twoway ///
(rarea u7a d7a  Years4,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u7 d7  Years4,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b7 Years4, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero4 Years4, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ebp1, replace

twoway ///
(rarea u8a d8a Years4,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u8 d8  Years4,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b8 Years4, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero4 Years4, lcolor(black)), legend(off) ///
title("Tight Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ebp2, replace

gr combine fig_ebp1 fig_ebp2, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("Ten-Year Term Premium", size(small)) ///

gr rename fig_baa, replace



* Combine all 4 IRFs
gr combine fig_ebp fig_cp fig_aaa fig_baa, ///
rows(2) cols(2) graphregion(color(white)) plotregion(color(white)) ///

********************************************************************************